function [] = gendata(hyp_mean, n)
% a script to generate data by a Gaussian process directly
% by Mark Norrish, 2011
% hyp_mean = [0;0];
if nargin < 2
  n = 125;
end
n_class = 3;

D = 2; % dimensionality
Data = rand(n, D);

func = @covSEard;

if nargin == 0
  hyp_mean = zeros(eval(func()),1);
end

%hyp_mean = [0 0 0]';
%hypmean=[0, 1:x,0]
%hyp_var = [1 1]';
%hyps = normrnd(repmat(hyp_mean, 1, n_class), repmat(hyp_var, 1, n_class)); % randomly generate hyps
hyps = -1.6*log(n_class)*ones(eval(func()), n_class); % a decent way of getting decentish data
hyps(:,1) = hyp_mean;
% 

sigma = 1e-4;
K = zeros(n,n,n_class);

for i = 1:n_class
  K(:,:,i) = func(hyps(:,i), Data, Data) + sigma*eye(n);
end

cholk = zeros(n*n_class);
for c = 1:n_class
  cholk(1+(c-1)*n:c*n,1+(c-1)*n:c*n) = chol(K(:,:,c));
end

Y = [];
while ~all(ismember(1:n_class, Y))

  F = cholk'*normrnd(zeros(n*n_class,1), ones(n*n_class,1)); % because mu is zero
  Fdash = reshape(F, n, n_class)';
  Y = (Fdash==repmat(max(Fdash),n_class,1));
  Y = sum(Y.*repmat((1:n_class)',1,n));

end
dlmwrite('gp_data', [Data Y']);
dlmwrite('gp_function', F);
dlmwrite('gp_hyps', hyps);
